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ABSTRACT 

We present results from the first three-dimensional radiation hydrodynamical calculations to 
follow the collapse of a molecular cloud core beyond the formation of the stellar core. We 
find the energy released by the formation of the stellar core, within the optically-thick first 
hydrostatic core, is comparable to the binding energy of the disc-like first core. This heats the 
inner regions of the disc, drives a shock wave through the disc, dramatically decreases the 
accretion rate on to the stellar core, and launches a temporary bipolar outflow perpendicular 
to the rotation axis that travels in excess of 50 AU into the infalling envelope. 

This outburst may assist the young protostar in launching a conventional magnetic jet. 
Furthermore, if these events are cyclic, they may provide a mechanism for intense bursts of 
accretion separated by long periods of relatively quiescent accretion which can potentially 
solve both the protostellar luminosity problem and the apparent age spread of stars in young 
clusters. Such outbursts may also provide a formation mechanism for the chondrules found in 
meteorites, with the outflow transporting them to large distances in the circumstellar disc. 

Key words: accretion, accretion discs - hydrodynamics - radiative transfer - stars: formation 
- stars: low-mass, brown dwarfs - stars: winds, outflows. 



1 INTRODUCTION 



More than four decades ago, |Larson| \ 1 969} performed the first nu- 
merical calculations of the collapse of a molecular cloud core to 
stellar core formation and beyond. These one-dimensional radia- 
tion hydrodynamical calculations revealed the main stages of proto- 
star formation: an almost isothermal collapse until the inner regions 
become optically thick, the almost adiabatic formation of the first 
hydrostatic core (typical radius ^ 5 AU and initial mass ^ 5 Mj), 
the growth of this core as it accreted from the infalling envelope, 
the second collapse within this core triggered by the dissociation of 
molecular hydrogen, the formation of the stellar core (initial radius 
^ 2 Rq and mass ^1.5 Mj), and, lastly, the long accretion phase 
of the stellar core to its final mass. Subsequent one-dimensional 
(e.g. Masunaga & Inutsuka 2000) and two-dimensional (Tschar- 
|nuter^l987 . Tscharnuter et al.,2009^ calculations have not changed 
this qualitative picture substantially, although the latter have al- 
lowed the disc-like structure of rotating first cores to be studied. 

The first three-dimensional hydrodynamical calculations to 
follow the collapse to stellar core formation were performed more 
than a decade ago by |Batej(^1998j). These calculations of rotating 
molecular cloud cores examined the non-axisymmetric evolution 
of the first core and the second collapse phase. If the first core 
was rotating rapidly enough it was found to be dynamically unsta- 
ble to a bar-mode leading to the formation of trailing spiral arms. 
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Gravitational torques removed angular momentum and rotational 
support from the inner regions of the first core, quickening the on- 
set of the second collapse and preventing fragmentation during the 
second collapse phase to form close binaries. Several subsequent 
studies have investigated this phenomenon in more det ail (jSaigo 



& Tom isaka'2006 , Saigo, Tomisaka & Matsumoto|2008[pachida^ 
Inutsuka & Matsumoto 2010)), some also including magnetic fields 
and finding outflows (Machida et al. 2005] |2006| ). However, all 
these calculations used barotropic equations of state rather than 
solve the radiation hydrodynamical problem. 

The first three-dimensional calculations including radiative 
transfer that followed collapse to the point of stellar core forma- 
tion (but not beyond) were IWhitehouse & BatelpOPg]), using the 
flux-limited diffusion approximation, and' Stamatellos et al.| ( |2007| ), 
using a radiative cooling approximation. Most recently, radiation 
magnetohydrodynamical calculations of cloud collapse have been 
performed ( [Tomida et al.|2010) , but were stopped before the onset 
of the second collapse. 

In this paper, we report results from the first three-dimensional 
radiation hydrodynamical calculations to follow the collapse of ro- 
tating molecular cloud cores beyond the formation of the stellar 
core. We find that the use of radiation hydrodynamics rather than a 
barotropic equation of state has little effect up until the formation 
of the stellar core. However, with radiative transfer, the energy re- 
leased by the formation of the stellar core has a dramatic effect on 
the surrounding disc and envelope and launches a temporary out- 
flow even in the absence of a magnetic field. 
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Figure 1. The time evolution of the maximum (central) density (left two panels) and temperature (right panel). For the left and right panels, the lines are for 
cloud cores with /5 = 0, 5 x 10~^, 0.001, 0.005, 0.01 from left to right. For the central panel, for each calculation the time in years has been zeroed when 
the density first exceeded 10"*^ g cm"*^ and the /3 = calculation is plotted with a solid line, ^ = b x lO"'^ with a dotted line, and the other calculations 
are indistinguishable. The free-fall time of the initial cloud core, fff = 1.8 x 10-*^^ s (56,500 yrs). Each calculation was performed with 10^ SPH particles. 



2 COMPUTATIONAL METHOD 

The calculations presented here were performed using a three- 
dimensional smoothed particle hydrodynamics (SPH) code based 
on the original version of |Benz| |l990 Benz et al. 1990), but sub- 
stantially modified as descri bed in ,3 ate et al.,^1995^),,Whitehouse, 
Bate & MonaghanI poos) ), [Whitehouse & Bate] ( |2006| ), |Price~ 
Bate p007X and parallelised using both OpenMP and MPL 

Gravitational forces between particles and a particle's near- 
est neighbours are calculated using a binary tree. The smoothing 
lengths of particles are variable in time and space, set iteratively 
such that the smoothing length of each particle h — 1.2(m/p)^/^ 
where m and p are the SPH particle's mass and density, respectively 
(see [Price & Monaghan 2007 for further details). The SPH equa- 
tions are integrated using a second-order Runge-Kutta-Fehlberg in- 
tegrator with individual time steps for each particle (IB ate et al. 
1995|). To reduce numerical shear viscosity, we use the M orris & 
Monaghan| ( |1997| artificial viscosity with varying between 0.1 
and 1 while fi^ — 2a^ (see also [Price & Monaghan|2005) ). 



2.1 Equation of state and radiative transfer 

We use an ideal gas equation of state p — pTlZ/ p, where T is the 
gas temperature, IZ is the gas constant, and p is the mean molecu- 
lar mass. The equation of state takes into account the translational, 
rotational, and vibrational degrees of freedom of molecular hydro- 
gen (assuming a 3:1 mix of ortho- and para-hydrogen; see |Boley] 
|et al.|2007] ). It also includes molecular hydrogen dissociation, and 
the ionisations of hydrogen and helium. The hydrogen and helium 
mass fractions are X = 0.70 and Y = 0.28, respectively. The con- 
tribution of metals to the equation of state is neglected. Two tem- 
perature (gas and radiation) radiative transfer in the flux-limited dif- 
fusion approximation is implemented as described by Whit ehouse] 
[eta l. (2005) and Whitehouse & Bate (2006), except that the stan- 
dard explicit SPH contributions to the gas energy equation due to 
the work and artificial viscosity are used when solving the (semi- 
)implicit energy equations to provide better energy conservation. 
We assume solar metallicity gas, using the interstellar grain opacity 
tables of Pollack et al. (1985) and the gas opacity tables of|Alexan-| 
|der| ( ^1975j (the IVa King model) (see , Whitehouse & Bate|2006| ). 



2.2 Initial conditions 

The initial conditions for the calculations are identical to those of 
|Bate| ( |1998| ). We follow the collapse of initially uniform-density. 



uniform-rotating, molecular cloud cores of mass M — 1 and 
radius R — 7 x 10^^ cm. The ratios of the thermal and rotational 
energies to the magnitude of the gravitational potential energy are 
a = 0.54 and /3, respectively. We have performed calculations with 
/3 = (not rotating), and /3 = 5 x 10"^, 001, 0.005 and 0.01. 

To satisfy the resolution criterion of [Bate & Burkertj ( |1997| l 
that the minimum Jeans mass during the calculation contains at 
least ^ 2A^neigh = 100 particles, we require at least 1 x 10^ equal- 
mass particles. To test for convergence, we performed calculations 
using 1 X 10^, 3 X 10^ 1 x 10^, and 3 x 10^ equal-mass SPH par- 
ticles. The calculations were performed on the University of Exeter 
Supercomputer, an SGI Altix ICE 8200. 



3 RESULTS 

The collapse of the molecular cloud cores up until the forma- 
tion of the stellar core proceeds in a manner very similar to 
that reported from previous three-dimensional calculations using 



barotropic equations of state (Bate 


1998 


JSaigo & Tomisaka|2006| 


ISaigo et al.||2008| |Machida et al. 


2010 


1. Fig. [1] gives the evolu- 



tion of the maximum density and temperature for the 10 par- 
ticle calculations with different initial rotation rates. The initial 
collapse is almost isothermal until the maximum density exceeds 



10- 



g cm .As the central regions of the cloud become op- 



tically thick, the gas begins to trap the radiation and the collapse 
enters an almost adiabatic phase where the temperature rises as 
the gas is compressed. This leads to the formation of a pressure- 
supported 'first core' ( |Larson|1969| ), which can be seen in Fig. [T] 
when the initial collapse stalls with central (maximum) densities 
^ 10~^^ g cm~^ and temperatures of ^ 70 — 120 K, depend- 
ing on the degree of rotational support (i.e. cores that rotate more 
quickly have lower central temperatures). Without rotation, the first 
core has an initial mass of 5 Jupiter masses (Mj) and a radius of 
^ 5 AU. However, with higher initial rotation rates of the molec- 
ular cloud core, the first cores become progressively more oblate. 
For example, with (3 — 0.005 before the onset of dynamical insta- 
bility, the first core has a radius of 20 AU and a major to minor 
axis ratio of ^A\\ (left panels, Fig.[2|. Thus, for the higher rotation 
rates, the first core is actually a disc, but without a central object. 
As pointed out by|Bate|([T998| and Machida et al.|p010 ), in these 
cases the disc actually forms before the star. 

The subsequent evolution of the first core depends on its rota- 
tion rate. Non-rotating and slowly rotating cores evolve as they ac- 
crete mass from the surrounding inf ailing envelope with their cen- 
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Figure 2. Column density images of the evolution of the first core, shock wave, and outflow for the /3 = 0.005 case (using 3 x 10^ SPH particles). The rapidly 
rotating first core is highly oblate (left panels) and undergoes a dynamical bar-mode instability (centre-left panels) which transports angular momentum out of 
the centre, transforms the core into a disc, and triggers the formation of the stellar core. The energy released by stellar core formation heats the inner regions of 
the disc and drives a shock wave through the disc (upper right two panels) and a bipolar outflow perpendicular to the disc (lower right two panels). Animations 
can be found at http ://www. astro .ex. ac .uk/people/mbate/Animations/Stellar/| . 



tral densities and temperatures increasing (Fig.[T] calculations with 
13 ^ 0.001). When the central temperatures of a first core exceeds 
^ 2000 K, molecular hydrogen begins to dissociate, absorbing en- 
ergy and leading to a second hydrodynamic collapse deep within 
the first core. The collapse continues until the molecular hydrogen 
has been complete dissociated at the centre of the core where upon 
a second pressure-supported core begins to form, the 'stellar' core 
([Larson 1969 ). The formation of the stellar core occurs just a few 
years after the onset of the second collapse, during which the max- 
imum density increases from 10~^ to ^ 0.1 g cm~^ and the 
maximum temperature increases from 2000 to > 60, 000 K. The 
stellar core is formed with a mass of Msc ~ 1.5 Mj and a radius 
of Rsc ~ 2 R© . Without rotation, the stellar core accretes the rem- 
nant of the first core in which it is embedded 10 yr and then 
accretes the envelope ( |Larson|1969| ). 

If the first core is rotating rapidly enough that its value of 
13 > 0.274, the core is dynamically unstable to the growth of non- 
axisymmetric structure ( |Bate|1998) ). For the particular initial con- 
ditions used here, this occurs for the (3 = 0.005 (see Fig.|2} and 
13 — 0.01 cases. The first core develops a bar-mode, and then spiral 
structure as the ends of the bar wind up (e.g. |Durisen et al.|1986{ 
|Bate|1998] ). The spiral structure removes angular momentum from 
the inner parts of the first core, the effect of which can be seen in the 
evolution of density and temperature in Fig.[T] Specifically, for the 
P = 0.005 case, the slow increase in central density and tempera- 
ture accelerate with the onset of the spiral structure at t = 1.05 tff. 
This substantially accelerates the evolution of the first core towards 
the second collapse, which would have taken much longer to reach 
without the angular momentum redistribution. It also inhibits frag- 
mentation during the second collapse phase because the gas un- 
dergoing the second collapse has less angular momentum than it 
would otherwise (Bate 1998 ). When the stellar core forms it is al- 
ready surrounded by a large disc (radius ^ 50 AU; Fig. [2]) which 



is the remnant of the gravitationally-unstable first core ( |Bate|1998l 
IMachida et al.|2010| ). 

For even higher rotation rates of the initial molecular cloud 
core (/3 ^ 0.02), the first core is actually a ring-like structure (e.g. 
[Norman & Wilson] 1 978 |[Cha & Whitworth|2003| ) which is prone to 
dynamical fragmentation into two or more objects, each of which 
evolves in a similar manner to the more slowly rotating first cores. 

As mentioned at the start of this section, all of this evolution 
with radiation hydrodynamics is almost identical to the evolution 
that is obtained when using a barotropic equation of state. The 
small differences that are found will be discussed in a separate pa- 
per. The main purpose of this letter is to discuss the evolution sub- 
sequent to the formation of the stellar core, which is qualitatively 
different to that obtained with a barotropic equation of state. 

3.1 The effect of stellar core formation 

When using a barotropic equation of state, the formation of the 
stellar core deep within the optically-thick disc has no effect on the 
temperature of the majority of the gas in the disc because the tem- 
perature is set purely according to the density of the gas. However, 
with radiation hydrodynamics, the situation is completely different. 

When the second collapse occurs and produces the stel- 
lar core, the gravitational potential energy that is released is ^ 
GMl/R,c = 4 X 10^^ erg. Since the stellar core is in virial equi- 
librium, approximately half of this energy is to be radiated away. 
Moreover, the stellar core rapidly begins to accrete from the first 
core, reaching a mass of 6 Mj in only a few years (i.e. more 
quickly than the dynamical timescale of the disc; see below). This 
increases the total energy released to 3 x 10^^ erg. However, 
the binding energy of the disc in the /3 = 0.005 calculation, for 
example, calculated just before the onset of the second collapse 
(i.e. the binding energy of the gas surrounding the stellar core that 
makes up the disc) is only 4 x 10^^ erg (which can be estimated as 
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Figure 3. A cross section down the rotation axis {x — z plane) of 
the (3 = 0.005 calculation with 3 x 10^ particles showing the den- 
sity (colour scale) and velocities (vectors). The bipolar outflow is clearly 
visible and has travelled 50 AU into the surrounding inf ailing enve- 
lope only 35 yr after stellar core formation. The shock wave travel- 
ling along the midplane of the disc has propagated 17 AU out into 
the disc in the same period of time. Animations can be found at 
|http://www.astro.ex.ac.uk/people/mbate/Animations/Stellar/| . 
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Figure 4. The mass of the stellar core (gas with density > 10~^ g cm~^) 
measured from the time of its formation for the /3 = 0.005 molecular cloud 
core performed using 10^ (long-dashed line), 3 x 10^ (short-dashed line), 
10^ (dotted line), and 3 x 10^ (solid line) particles. With low resolution the 
radiative feedback stops accretion onto the stellar core, but for the highest 
resolution the accretion continues at a low level (10~^ Mq yr~^). 



- GMi/Rd with Md ^ 0.18 M© and taking a mean 'spherical' 
radius of Rd ~ 15 AU). Thus, the disc suddenly finds itself irra- 
diated from the inside by an energy source emitting a substantial 
fraction of the binding energy of the disc itself. Because the disc is 
extremely optically thick, this energy is temporarily trapped in the 
centre of the disc and heats the gas dramatically, sending a weak 
shock wave (Mach number ^ 1) out along the midplane of the disc 
(initial speed ^ 2 km s~^). However, perpendicular to the disc, the 
effect is even more dramatic. Because there is less material along 
the rotation axis, the hot gas finds it easiest to break out in this di- 
rection and a bipolar outflow is launched. Whereas the wave within 
the disc decays as it travels leaving the bulk of the disc gravitation- 
ally bound, the gas forming the bipolar outflow has velocities in 
excess of 8 km s~^ and travels out into the infalling envelope to 
distances in excess of 50 AU in less than 50 years (Figs.[2]and|3]). 



The energy released by stellar core formation is somewhat de- 
pendent on the resolution of the calculations because the heating 
and outflow inhibit accretion onto the stellar core, which in turn 
results in less energy being released. Thus, the higher the resolu- 
tion, the quicker the radiative feedback is able to inhibit the accre- 
tion and, therefore, the lower the initial mass of the stellar core. 
In Fig. [4] we plot the mass with density > 10~^ g cm~^ versus 
time for the /3 = 0.005 case with resolutions ranging from 10^ to 
3x10^ SPHparticles.lt can be seen that the mass of the stellar core 
at which the feedback dramatically curtails the accretion decreases 
from ^ 11 to 6 Mj with increased resolution. For a few years, the 
stellar core grows at a rate of ^ 10~^ M© yr~^ (even with the high- 
est resolution). With low-resolution (^ 1 x 10^ SPH particles), the 
accretion onto the stellar core actually ceases entirely during the 
launching of the outflow. However, using 3 x 10^ SPH particles, 
we find that the stellar core does continue to grow during this time, 
but at a much reduced accretion rate of 1 x 10~^ M© yr~^ (mea- 
sured between 20 and 50 years after stellar core formation; Fig.|4]). 
Although convergence with increasing resolution is relatively slow, 
a strong outflow is launched in all cases. 



4 DISCUSSION 

These first calculations raise many questions to be followed up in 
the future. First, although the outflow seen in the calculations dis- 
cussed here is transient and cannot be the origin of proto stellar jets 
or outflows, such heating upon stellar core formation might aid in 
the launching of a magnetic jet by helping to clear material perpen- 
dicular to the disc and combining substantial thermal support with 
magnetic support. Therefore, using future radiation magnetohydro- 
dynamical calculations it will be important to assess the relative 
roles of the radiative impact of stellar core formation and magnetic 
fields in the launching of the proto stellar jet. Similarly, it would be 
desirable to resolve the accretion shock at the surface of the stellar 
core and to use a radiative transfer method that treats the optically- 
thin regime better than the flux-limited diffusion approximation. 

Second, although the dramatic decrease in the stellar core ac- 
cretion rate is also a transient, the question arises as to whether 
such a situation might be cyclic. One can imagine that after the 
disc settles down again, the accretion rate on to the stellar core 
may increase, leading to a repeat of the rapid accretion, the genera- 
tion of a large accretion luminosity, the strengthening of an outflow 
and heating of the disc, and the consequential decrease in the stel- 
lar core's accretion rate. The driving mechanism for the increased 
stellar accretion rate could be the redevelopment of spiral density 
waves in the disc as it recovers from the previous shock wave 
(Fig. [2]). Such gravitational instabilities have been invoked to ar- 
gue that accretion in the early Class and I stages of star formation 
may be highly variable ( ,Vorobyov & Basu|2005||200 6 ). The stellar 
accretion luminosity may then provide a way to temporarily switch 
the gravitational instability off. Another possible mechanism for 
variable accretion might be a mismatch between the accretion rates 
given by gravitational instability at large radii and the magnetoro- 
tational instability at small radii ( [Armitage et al.|200T| |Zhu et al^ 
[2009). If most of a star's mass is accumulated in short bursts of 
intense accretion, this can potentially explain the observation that 
protostars are observed to be less luminous than expected for steady 
accretion (see |Kenyon et al.|1990[|Hartmann & Kenyon|1996| ), and 
also the apparent age spread of young stars in the Hertzsprung- 
Russell diagram (see |Baraffe, Chabrier & Gallardo|2009] ). 

Finally, there is the age old problem of how the chondrules 
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found in meteorites were formed ( [Grossman et al.|p"988| ) as their 
production is thought to require high temperatures (> 1500 K). 
Although many theories have been advanced in order to explain 
chondrule formation (see, for example, the review by Boss 1996 ) , 
there is still no widely accepted mechanism. In this context, it is of 
interest to note that the temperatures within ^ 3 AU of the stellar 
core at the launching of the outflows exceed 1500 K. While most of 
this material is no doubt be accreted by the star, a fraction produces 
the outflows which reach beyond 50 AU. If chondrules were formed 
in this outburst event, some might be entrained in the outflow and 
eventually fall back into the circumstellar disc at large distances 
from the star. Again, if such events were cyclic this would help be- 
cause chondrule formation is thought to have occurred over several 
million years (Scott,2007j ). 



5 CONCLUSIONS 

We have presented results from the first three-dimensional radiation 
hydrodynamical calculations to follow the collapse of a molecular 
cloud core beyond the formation of the stellar core. We find the 
evolution before the formation of the stellar core is very similar to 
that found in the past using barotropic equations of state. In partic- 
ular, the evolution of the first hydrostatic core which, depending on 
its rotation rate, may be dynamically unstable to the growth of non- 
axisymmetric perturbations and the generation of spiral structure, 
is very similar to that found in barotropic calculations. As found 
in earlier calculations, a rapidly rotating first core actually evolves 
into a disc so that the disc actually forms before the stellar core. 

However, we find that the evolution following the formation of 
the stellar core is qualitatively different in radiation hydrodynam- 
ical calculations. In barotropic calculations, the formation of the 
stellar core deep inside the first core (or disc) has no effect on the 
surrounding disc because the temperature of the gas is simply set 
by the density of the gas. However, with radiation hydrodynamics, 
the energy released by the formation of the stellar core within the 
optically thick disc is similar to the binding energy of the disc. This 
heats the inner regions of the disc, drives a shock wave outwards 
through the disc, dramatically decreases the accretion rate on to the 
stellar core, and launches a bipolar outflow perpendicular to the ro- 
tation axis that can travel in excess of 50 AU out into the infalling 
envelope in less than 50 years. 

We speculate that such outflows may assist the young proto- 
star in launching a conventional magnetic jet by clearing a path 
perpendicular to the disc and adding substantial thermal pressure 
to the force provided by the magnetic field. It may also be that such 
events are cyclic, occurring every time the accretion rate onto the 
protostar exceeds a certain level rather than simply being a one-off 
event associated with the formation of the stellar core. If so, such 
events may provide the mechanism for intense bursts of accretion 
separated by long periods of relatively quiescent accretion which 
may be necessary to solve the protostellar luminosity problem and 
the apparent age spread of young stars. Finally, such outbursts may 
provide another mechanism for the formation of chondrules in me- 
teorites and the associated outflows may be able to transport them 
to large distances in the circumstellar disc. 
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